Momentum distribution of the uniform electron gas: 
improved parametrization and exact limits of the cumulant expansion 
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The momentum distribution of the unpolarized uniform electron gas in its Fermi-liquid regime, 
n(k, r s ), with the momenta k measured in units of the Fermi wave number fcp and with the density 
parameter r 3 , is constructed with the help of the convex Kulik function G(x). It is assumed that 
n(0,r s ), n(l ,r s ), the on-top pair density g(0,r s ) and the kinetic energy t(r s ) are known (respec- 
tively, from accurate calculations for r s = 1, ...,5, from the solution of the Overhauser model, and 
from Quantum Monte Carlo calculations via the virial theorem). Information from the high- and 
the low-density limit, corresponding to the random-phase approximation and to the Wigner crys- 
tal limit, is used. The result is an accurate parametrization of n(k,r s ), which fulfills most of the 
known exact constraints. It is in agreement with the effective-potential calculations of Takada and 
Yasuhara [Phys. Rev. B 44, 7879 (1991)], is compatible with Quantum Monte Carlo data, and is 
valid in the density range r s < 12. The corresponding cumulant expansions of the pair density and 
of the static structure factor are discussed, and some exact limits are derived. 

PACS numbers: 71.10.Ca, 05.30.Fk 



I. INTRODUCTION 

In solid state theoryS and quantum chemistry,i the 
phenomenon of electron correlation and some of its de- 
tails are hidden Jja-tLie reduced densities and reduced 
density matricesBtatffl and their cumulants.DEl For the 
ground state of the uniform electron gas (jellium) of den- 
sity p = 3/(47rrj) (in a.u.), these quantities are the pair 
density g(x, r s ) and the momentum distribution n(fc, r s ), 
where k is measured in units of the Fermi wave num- 
ber k-p = (S^p) 1 / 3 and x is the scaled interelectronic 
distance, x = k~pr\2- 

Besides its relevance in the understanding of many ef- 
fects in simple metals and semiconductors, the jellium 
model plays a crucial role in providing input quanti- 
ties for approximate approaches to the many-electron 
problem of nonuniform density. Different approximate 
schemes, in fact, often need different quantities from jel- 
lium. As an example, density functional theory (DFT) 
uses the exchange-correlation energy of the uniform elec- 
tron gas for the widely used local density approximation 
(LDA). For building nonempirical beyond-LDA function- 
als for use in DFT, the pair density of jellium is often 
needed. El 

In such applications, the quantities from the uniform 
electron gas must be available in the form of analytic 
expressions. Since the jellium model is not exactly solv- 
able, analytic expressions are built by interpolating be- 
tween known exact limits, and by fitting the Quantum 
Monte Carlo (QMC) data, when available. Relevant ex- 
amples are (i) the correlation energy used ift-LDAi imple- 
mentation, built by using functional formscJ : l!3tij which 
includes exact limite-jand interpolate the QMC data of 
Ceperley and Alder;!! 3 ! (ii) the pair density, built by com- 



bining exact properties and fitting to QMC data 



or by simply interpolating between exact limits;E3 (i 
the static local field factor, parametrized by fittinL 
QMC data of the static response;ta (iv) the dynami 
cal local field, factors, built by using many known exact 
constraints O All these parametrized quantities are not 
strictly "exact", but are considered to be the closest to 
the true quantities, in the sense that this is the best one 
can presently obtain, and in the sense that they are ac- 
curate enough for the purpose for which they are needed. 

In recent years, there has been increasing interest in 
a particular approach to the many-electron problem of 
nonuniform density, the so called density matrix func- 



ses the one-body density 
1 J2 ' Building a "local ap- 



tional theory (DMFT), wi 
matrix as basic variable.^ 
proximatkm" for DMFT is not an easy task: in a first 
attempt,! 2 ^ the momentum distribution n(k, r s ) of jellium 
has been used as input. Besides this important appli- 
cation, for which a reliable parametrization of n(k,r s ) 
is needed, there are other reasons to focus on the mo- 
mentum distribution of jellium. The definitions of "ex- 
change" and "correlation" in DMFT are different from 
those in DFT and in Hartcee-Fock-like approaches: in 
DMFT the cumulant partLra of the pair density rises 
to be a key quantity. An accurate parametrization of 
n(k,r s ) at metallic densities allows to extract the cumu- 
lant pair de.ns.ity.-Qf jellium, since the whole pair density 
is available otJ'EJ The cumulant pair density can then 
be compared with recent -attempts to calculate it in a 
high-density electron gasJ^'LJ and can be diagonalizcd 
in terms of "cumulant geminals^i^ftalog of "Overhauser 
geminals" for the pair densityE§EZI£3) . Also, the study 
of exact limiting behaviors of the cumulant pair density 
is of great interest, since some of these limits could be 



2 



approximately valid in nonuniform systems. 

The momentum distribution of the uniform electron 
gas is also useful for the calculation of the exchange- 
correlation correction to Compton profiles computed 
within LDA-DFTB3 Note also that n(k,r s ) determines 
the part of the local-field factor beyond the random- 
phase approximation (RPA) which takes into account the 
change in occupation numbers in the Lindhard function 
[seeEq. (29) of Ref.[j]J. 

The qualitative behavior of n(k, r s ) is the following. It 
starts at k = with a value n(0, r s ) < 1 and decreases 
with increasing k. For fc < 1 it is concave. Then in 
the Fermi-liquid regime at k = 1 there is a finite jump 
(Fermi gap) from n(l~,r s ) to a lower value n(l + ,r s ) = 
n(l~,r s ) — z-p(r s ) with logarithmic slopes at both sides 
of k = 1. For k > 1 (correlation tail) n(k,r s ) is convex, 
and vanishes for k — > oo. For r s = (ideal Fermi gas) 
it is n(k,0) = 9(1 — k), where 8(x) is the Heaviside step 
function. Thus, the quasiparticle weight zp(r s ) starts 
with zf(0) — 1 and decreases with increasing interaction 
strength r s . At large r s the electrons form a Wigner 
crystal with a smooth n(k, r s ). r s <C 1 and r s 3> 1 
are the weak- and strong-correlation limits, respectively. 
For intermediate values of r s a non-Fermi liquid regime 
may exist with z-p = 0. In such case n(k,r s ) would be 
continuous vs. k, with a nonanalytical behavior at k = 1. 

In Rcf. |30[ the idea of using the convex Kulik func- 
tion G(x) to parametrize the two branches (k < 1 and 
k > 1) of n(fc,r s L-is sketched. The function G(x) ap- 
peared in Kulik'sEil RPA analysis of n(fc,r s ) near the 
Fermi edge (|1 - k\ « 1), see Eq. (|A4|). G(x) be- 
haves as Co + cixlnx for small x (see Appendix and 
Fig. [I]) , which corresponds to the correct nonanalytic be- 
havior of n(k,r s ) near the Fermi surface. So, suppos- 
ing that the value at the centre, n(0, r s ), and the values 
at the Fermi edge, n(l~,r s ) and n(l + ,r s ), are known, 
it should be possible to represent n(k,r s ) in terms of 
G(x), with suitable prefactors and with suitable scaling 
(squeezing and stretching) of its argument. In this way, 
n(k,r s ) becomes a functional of no(r s ) = n(0, r s ), and 
of n±(r.) - n(l± r s ), and can be designed to yield the velopments are reported in the last Sec. |V3 

proper normalization and the correct kinetic energy t(r s ), 
which follows from the-total energy e(r s ) = t(r s ) + v(r s ) 
via the virial theoremE3 In R^frj30[ the input data from 
Takada and Yasuhara (TY)BE5for n (r s ), n±(r s ) and 
t(r s ) at r s = 1, 5 have been used, together with the 
on-top pair density go(r s ) = g(Q,r s ) [which determines 
the large-fc behavior of n(k)] from Ref. ^?]. The result is 
a field n(k,r s ) for r s € [1,6] which is correctly concave 
for k < 1, convex for k > 1, and with a Fermi gap zf(t s ) 
at k = 1 decreasing with growing r s . The attempt to 
extend this procedure for r s 6 [6, 10] failed: n(k < l,r s ) 
is no longer concave for r s > 6. 

Here, an improved version of the parametrization of 
n(k, r s ) in terms of the Kulik function G(x) is presented. 
Our parametrized momentum rdjjstpbution is in good 
agreement with the TY values,! 3 ti and is valid in the 
range of densities r s < 12. Previous parametrizations of 



n(k, r s )SS0@ used the Quantum Monte Carlo (QMC) 
data of Ref. |3^ as an input. However, QMC data are 
presently only available for 0.4 < k < 0.9 and 1.1 < k < 
1.5, thus not providing information about n(k,r s ) near 
the centre, k — 0, and at the Fermi edge, k = 1. In 
these last regions, in fact, different parametrizations of 
the same_.QMC data can be rather different from each 
other .EEIeSlZI Our construction of n(k, r s ) uses informa- 
tion from the jefEeiptive-potential calculations of Takada 
and Yasuhara,E3o from the high- and low-density limits 
of n(k,r s ), corresponding to RPA and the Wigner crys- 
tal (|WC) limit, and iom accurate parametrizations of 
t(r s )E3 and of go (r s )E3 In the regions where QMC data 
are available, our n(k, r s ) is compaJihIe |[ J5z | ith them. Also, 
with respect to previous works j23-E3'EZrE3 the functional 
form used here satisfies more exact limits. In particu- 
lar, the logarithmic behavior jat the Fermi edge is taken 
into account for the first time-Ea Notice that it causes the 
logarithmic divergence of t(r s — ► 0)o 

Using our n(k,r s ), the moments (k v ) : the correla- 
tion entropy, and the one-body reduced density matrix 
f(x,r s ) are evaluated. The latter appears in the cumu- 
lant partitioning of the pair density g(x,r s ). The static 
structure factor S(k,r s ), related to the pair density via 
Fourier transform, the particle-number fluctuations in 
fragments AA r n(r s ), and the potential energy v(r s ) are 
discussed in terms of their cumulant partitioning, and 
some exact limits are derived. Finally, by meanSj-af an 
accurate model for the spin-resolved pair density,E£l the 
cumulant part of g(x, r s ) is extracted. 

The paper is organized as follows. In Sec. the known 
sum rules and limiting cases for n(fe, r s ) are reported, and 
they are used in Sec. HI to build up our parametriza- 



tion of the momentum distribution via the Kulik func- 
Sec. 



tion. 



IV is devoted to the calculation and discussion 



of different measures of the correlation strength, of the 
one-matrix, and of the cumulant expansion of the pair- 
density. In Sec. we study the cumulant partitioning 
of the static structure factor, of the density fluctuations 
and of the potential energy. Conclusions and future de- 



ll. SUM RULES AND LIMITING CASES 

How is n(k, r s ) defined? Starting from the many-body 
wave function ^(l, N), the one-body reduced density 
matrix (1-matrix for short) results from the N — 1 con- 
traction 



7(1|1') = 



dl...dN 

M 



d2...dN 

(iv - Ty. 



*(1,2,...,A0**(1',2,...A0,(1) 



|*(l,...,iV)| 2 = l, 



with the notation 1 = (r^ai). For the uniform electron 
gas, 7(1|1') = p<S CTlCT ; /(AjfIzj -ri|,r a ) defines the dimen- 
sionless 1-matrix f(x,r s ). Then its Fourier transform is 
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FIG. 1: The Kulik functi on G (x) appearing in the RPA anal- 
ysis of n(k,r a ), see Eq. (A4). At the origin G(x) has a fi- 
nite value, G(0) = 3.35 3337 [Eq. (A7)], and a logarithmic- 
divergent slope, see Eq. 



(A6) 



the momentum distribution 

„3 



n(k,r s ) = — dx 
1 Jo 



3 s in kx 

-kx~ f{x ^ 



(2) 



where dx 3 = 3/4-Kd 3 x = 3x 2 dx, and a = (4/97T) 1 / 3 . 
n(k,r s ) can be calculated using,perturbation theory, di- 
rectly with Green's functions] 3 .!! see Figs, la and lb of 
Ref. E4, or via the Hellmann-Feynman theoremcil as the 
energy-derivative n(k,r s ) — SE/dek, supposed E is (per- 
turbatively) known as a functional of Sk = hk 2 /2m and 
v q = 47re 2 / q 2 .E 2 ] Perturbative methods only work for 
high-densities, r s < 1. At metallic and lower densities, 
other techniques pajust be used, namely, the effective- 
potential methodJ23 which combines perturbation the- 
ory (Green's functions) with the Fermi rlwn ernetted 
chain approach, and the QMC simulations .E§c3 A more 
complete list of references concerning calculations and 
parametrizations can be found in Ref. |30|. 

n(k,r s ) has to satisfy the condition < n(k,r s ) < 
1 (which guarantees the ensemble A^-representability of 
the 1 -matrix), and the sum rules (k in units of kp, and 
energies in ryd) 



dk 3 n(k,r s ) = 1, (3) 



(«0 2 Jo 



dk 3 n(k,r s )k 2 = t(r s ), (4) 



where t(r s ) can be written as the sum of the kinetic en- 
ergy of the free Fermi gas, 3/5 (ar s ) -2 , and of the kinetic 
energy of correlation, t COIT (r s ). For r s < 1, t COTI (r s ) is 
known from RPA ajpiLicpm the lowest-order exchange 
diagram beyond it;Cj'OC3 for a summary see Eq. (3.25) 
and Figs, la and lb of Ref. ^4j At larger r s , t COIT (r s ) can 
be obtained via the virial ptheoremE3 from parametrized 
QMC correlation energies.EBl 

The large-fc behavior of n(k,r s ) is determined by the 
kinks in the many-body wave- function which occur when- 
ever two electrons are at contact or "on-top" (coalescing 



cusp properties) 



n(k —+ oo, r s ) 



C(r s 



k 8 \ k 10 



9?r 



2 {ar s ) 2 g (r s ), 



(5) 



where go(r s ) = g(0,r s ) is the on-top value of the pair- 
distribution function. In the r s — > QJiaait, go(r s ) can be 
obtained from perturbation theory,c2lE3 and at larger r s 
it has been calculated— bit, solving an effective two-body 
Schrodinger equation.EZHEl .— , 

In a normal Fermi liquid, the momeptum distribu- 
tion has a discontinuity and infinite slopesH at the Fermi 
edge, k = 1, 

n(fc^l~,r s ) = n-(r s ) - A(r a )(l - fc)ln(l -k) + 

0(1 -k), (6) 

n{k->l+,r s ) = n + (r s ) + A{r s )(k-l)\n(k-l) + 

O(k-l). (7) 

In the following, A(r s ) is referred to as the Fermi edge 
coefficient. In the small-r s limit, n±(r s ) and A(r s ) are 
known from RPA (see Appendix |X|). In the low-density 
or WC limit, the Fermi gap disappers [n_(r s ) — » n + (r s )], 
the infinite slopes at the Fermi edge may also vanish 
[A(r s ) -> 0]. 

Near the centr^ ^ — > 0, n(k,r s ) should behave 



quadratically, 

n(k -> 0, r s ) = n (r s ) + B(r s )k 2 + (9(fc 4 ). (8) 

A simple argument in favour of Eq. (||) is that it holds 
both in the high- and in the low-density limit (see Ap- 
pendix [A| and |b|) . 

When r s — >, Q . e xact results for n(k, r s ) are known by 
means of RPABEI (Appendix [Aj). In the RPA treatment, 
the Kulik function G(x) appears,EJ see Eqs. ( A3)-([A5|) 
and Fig. |l|. G(x) will be used in the next section to build 
up a parametrized n(k, r s ) which satifies Eqs. (§)-(§)■ 

In the low-density or strongly correlated limit, r s — > 
oo, the electron gas undergoes Wigner crystallization 
(see, e.g., Refs. |5^j5^ ). A simple model for the mo- 
mentum distribution in such regime is reported in Ap- 
pendix pl 



III. IMPROVED PARAMETRIZATION OF 

n(k,r a ) 

The momentum distribution in terms of the Kulik 
function G(x) of Fig. |l| is parametrized as follows. For 
k < 1 we use the ansatz 

n<(fc,r s ) = no - [ "° ~ ?~ ] G[ X< (k,r a )], (9) 



G(0) 



while for k > 1 we use 

n>(fc,r s ) = 



G(0) 



G[ x> (k,r.)], 



(10) 



4 



with x < (k,r s ) and x > (k,r s ) equal to 



x < (k,r s ) = a 



x > (k,r s ) = a 



ar s G(0) (1 - k) 
2w 2 [n - y/4ar s /iv 

vr 2 /7r(l-ln2) [n -n_] (1 - 
^Ty3 F"(0) G(0) a - 

ar s G(0) (k - 1) 



27r 2 7i + y/AarJn 

/37r(l-ln2) n+ tt 
' flo G(0) 4ar, 



(fc-1) 4 . (12) 



Here F"(0) = 17.968746 [see Appendix [A], Eq. Q], 
and the r s dependence of a, 6, no, n± and go is not 
explictly shown for shortness. These constructions are 
such that n<(fc) — > no,n^ for k — > 0, 1~, respectively, 
and n > {k) —* n+,0 for fc — » l + ,oo, respectively. The 
behavior of the Kulik function for small and large argu- 
ments (see Appendix 3) ensures the exact asymptotic 
expansion of Eqs. (|])-(|) near the centre, near the Fermi 
surface, and for large k. 

The parameter a(r s ) determines the Fermi edge coef- 
ficient A(r s ) of the |1 — k\ In 1 1 — k\ term at the Fermi 
surface, 



n(k — > 1^, r s ) 



„ ±(r ,) ± „w(^) 1,2 lQ + ^)x 

|l-fc|ln|l-fc| + 0(|l-fc|). (13) 



The parameter b(r s ) determines the curvature B(r s ) of 
Eq. (0) at the centre, k = 0, 



n(fc -> 0,r s ) = n (r s ) 



F"(0) 



b(r s 



■ 2 + 0(k 4 ). 

(14) 
0) = 1 and 



For small r s (RPA - Appendix [A|) it is a(r s 
b(r s -+ 0) = 1. 

In the preliminary version of Ref. |3^, another (but sim- 
ilar) ansatz was introduced, and it was chosen b(r s ) = 1. 
Two different functions, a<(r s ) and a>(r s ), for the co- 
efficient of |1 — As| In |1 — As | at k = 1~ and k = 1 + were 
fixed by the sum rules of Eqs. (0) and (Q) (with t CO rr(?*s) 
from Ref. [l2]). The values no^s) an( i rl ±( r s) were taken 
from the TY data (available for r s — 1, 5). The on-top 
value go(r s ) was taken from Ref. 27 

In our improved ansatz of Eqs. (|9])-(|l2|), we set 
a < (r s ) = a>(r s ) = a(r s ) (in agreement with Ref. |5^; 
also Figs. 7 and 8 of Ref. confirm this), and we use 
again t COIT (r s ) from Ref. ^2f and go(r s ) from Ref. |27]. 
Since we want to extend our results in the density range 
6 < r s < 10, where there are no data available for no(r s ) 
and n±(r s ), we extract information from the extreme 
low-density limit (Wigner crystal - see Appendix [^) , by 
following an oversimplified version of the idea presented 
in Ref. 

We first build no(r s ) by using a functional form which 
recovers the exact high-density limit, includes the Wigner 
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FIG. 2: Parametrized n(0, r s )Hjse|id line), compared to the 
Takada-Yasuhara (TY) valuesEJ-N The high-density or RPA 
limit and the Wigner crystal (WC) limit are also shown. 
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FIG. 3: Parametrized coefficient of the k term near the centre 
(* -» 0), B(r s ) = -sj£^2I [w^] 2 - The high-density or 
RPA result, b(r s 0) = 1, and the Wigner crystal (WC) 
limit are also shown. 



crystal behavior as r s — > oo, and has some free parame- 
ters to be fitted to the TY data. The result is reported 
in Fig. |2|, together with the high- and low-density curves. 
It is given by 



n (r a ) 



li+ 2,4- 5/2 



(15) 



with tt = 0.003438169, t 2 = 0.00725313666, t 3 = 
0.014900367, t 4 = 0.00113244364 (^ - t 3 agrees with the 
RPA value -(^) 2 4.1123 = -0.01146). We then build 
the parameter b(r s ) by a simple interpolation between 
the high- and low-density limit of the curvature at the 
centre (see Fig. 0). The result is 



b{r s 



1 + 0.0009376925 rf /4 ) 



1/2 



(16) 



Finally, the values at the Fermi edge, n±(r s ), and the 
coefficient of the infinite slope at the Fermi edge, a(r s ), 
are obtained by fitting the TY values for ?i±(r s )— while 
imposing the normalization and the kinetic energyll3 sum 
rules of Eqs. (0) and (Eh. The results are parametrized 
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FIG. 4: Upper panel: parametrized n(l ± ,r s ) (sqlijWines), 
compared to the Takada-Yasuhara (TY) valuesJ— I'l— I The 
high-density or RPA limit and the Wigner crystal (WC) 
limit are also shown. Lower panel: value of the Fermi gap 
zf(t s ) = n(l~ , r s ) — n(l + , r s ) as a function of r s ; (dip-present 
parametrization is compared with the TY results.Eja 



momentum distribution calculated with 
1,2, 10 (solid lines). In the Wigner 



FIG. 6: Upper panel 
Eqs. (§)-((l|) for r s ■. 
hmit, n(k,r s 3> 1) is calculated with Eq. (Bl) (dashed line, 
corresponding to r s = 75). Lower panel: ccp%p«rison of the 
present work with the Takada-Yasuhara (TY)OLj momentum 
distribution and the QMC calculation of Ref. ^ for r s = 5. 
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FIG. 5: Parametrized coefficient of the infinite slope at the 
Fermi edge, A(r s ) — a(r s ) (^f-) 1 ^ 2 0.63. The present result is 
compared with the RPA value, a(r s — > 0) = 1. Notice Ref. 



with the inclusion of the high- and low- density limits, 
and are equal to 



n-(r„) = 



1 + vi r s + v 2 r 2 + v 3 



1 + V4r s + v 5 rj +v e r^+ V7 r s 



15/4 ' 



(17) 



with vi = -0.0679793, v 2 = -0.00102846, v 3 
0.000189111, v 4 = 0.0205397, w 5 = -0.0086838, v 6 



6.87109-10- 



v 7 



4.868047 • 10" 5 {v x - v A agrees with 



the RPA value 



■3.3533 = -0.088519), and 



qi r s 



I , 1/2 i 7/4' 

l + q2r s +qsr s ' 

with q x = 0.088519 (from RPA), q 2 
0.022786335; 



(18) 



0.45, q 3 



a(r s ) 



1 



1/4 . 1/2 

Pir s +P2r s 



1 + P3 rl /4 + Pi rf " + P5 r s + p 6 r 



1/2 



(19) 



with pi = -78.8682, p 2 = -0.0989941, p 3 = -68.5997, 
p 4 = 38.1159, p 5 = -17.6829, p 6 = -0.01136759. Our 
parametrized n(k,r g ) breaks down at r s > 12 [in the 
density range 12 < r s < 16, n{k > l,r s ) is no longer 
convex, and for r s > 16 the unphysical result n_ < n + 
is obtained when the sum rules of Eqs. (§) and (Q) are 
imposed] . 

In the upper panel of Fig. ^, we show the functions 
n±{r a ), together with the TY values, and the high and 
low-density limits (here the r s — > oo limit is considered to 
be the in flex ion point of the WC momentum distribution, 
see Eq. (B4)). As said, our model is only valid for r s < 



12, so that n_(r s ) and n + (r s ) at densities lower than 
r s = 12 are no more obtained from the constraints of 
Eqs. (||) and (Q). Thus, the strange behavior of n_(r s ) 
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FIG. 7: The moments (k u ) of n(k,r s ) for r 3 = 3 and r s = 10. 
The corresponding results for the noninteracting gas (r s = 0) 
and for the Wigner crystal at r 3 = 75 are also reported. 



at r s ~ 16 does not affect our results. Also, the scheme 
presented here for the transition between the metallic 
and the extreme low-density region is oversimplified and 
must not be regarded as rigorous or reliable: we did not 
take into account the transition to the partially polarized 
electron gas (which affects the r s > 50 densitiesE^I) , as 
well as many other features. However, our results seem 
to be reliable in the relevant density range r s ^,12, and 
the simple picture of the upper panel of Fig. is only 
a "naive suggestion". In the lower panel of Fig. |J, we 
compare our parametrized z-p(r s ) with the TY results. 
In Fig. ||, we report the r s dependence of the Fermi edge 
coefficient A(r s ). Finally, in Fig. ^, we present in the 
upper panel our parametrized n(k,r s ) for 1 < r s < 10, 
and in the lower panel we compare our result with the 
TY n(k, r s ) and with the QMC data of Ref. || for r s = 5. 



IV. MOMENTS, CORRELATION ENTROPY, 
1-MATRIX, AND CUMULANT EXPANSION 



With the now available 
n(k,r s ), its moments 



momentum distribution 



dk 3 n(k,r s )k", 



(20) 



can be evaluated in addition to the normalization for v = 
and the kinetic energy for v = 2 [Eqs. (||) and @]. 
They are shown in Fig. |^ for r s = 3 and 10, together 
with the r s = (ideal Fermi gas) and with the WC result 
(r s = 75). The expression 



(Ai) 2 



1 



(ar s 



(21) 



(measured in ryd 2 ) describes the fluctuation of the kinetic 
energy. The moments (k 2 ) and (fc 4 ) determine the small- 
x behavior of the 1-matrix, see Eq. (p5|). 

the entropy- like expression s(r s ) = 



In Refs. 57 58 



FIG. 8: The particle-hole symmetric correlation entropy 
[Eq. (HU)], the normalization of the cumulant pair density 
[r.h.s. of Eq. Q34 )_]_, and the correlation tail normalization 

^corrtail 

(r s ) [Eq7"fl23|)] as a function of r s . 



has been used as a measure of the correlation strength.L£l 
Here the expression 



s p h(r s ) 



dk 3 (-l){n(fc,r s )lnn(fc,r s ) + 



[l-n(*c,r,)] ln[l-n(fc,r.)]} 



(22) 



is introduced as an alternative with the understanding 
that n(A;, r s ) and 1 — n(k, r s ) are the probabilities for the 
momentum state k to be occupied (with spin up and spin 
down) and empty, respectively. The entropy of this prob- 
ability 'distribution' is just the integrand of Eq. (p2|), and 
Sph{f s ) is the sum of all these entropies. Notice its invari- 
ance under the exchange n(k,r s ) <-> 1 — n(k, r s ), which 
is referred to as particle-hole symmetry in the Reduced- 
Density-Matrix community. This symmetry is an intrin- 
sic property p£ the correlation energy as a functional of 
the l-matrix.Eil s p h(r s ) is plotted in Fig. ||. Another mea- 
sure of the correlation strength is the correlation-tail nor- 
malization 



"corrtail = J dk 3 n(k,r s ), (23) 

also reported in Fig. M. For large r s , the Fermi edge 
disappears, z-p(r s ) = 0, and also any relict of it, A(r s ) = 
0, then the inflexion point of n(k,r s ) vs. k may serve in 
Eq. ( p3| ) as the lower limit. 

With n[k, r s ) also the (dimensionless) 1-matrix 

r°° oi n h T 

f(x,r s )= dk 3 - n(k,r s ), x = k F \r- r'| (24) 

Jo kx 

is available as the inverse of Eq. (||) . It has the small- a; 
behavior 



f(x< l,r s ) = 1 



(fc 2 ) , , (fc 4 ) 4 



-(lnn(A;, r s )) as a function of the interaction strength r s 



3! 
1 2 

5! 9^ 



5! 

(ar s ) 2 g (r s )-x 5 +0(x 6 )(25) 
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FIG. 9: 1-matrix f(x,r s ) for r s = (ideal Fermi gas), for 
r s — 10 (present model), and for r s = 75 [Wigner crystal 
(WC) limit, Eq. (ImJ)]. 



and the large- a; asymptotics (Friedel oscillations with re- 
duced amplitudes) 



/(x> l,r.) 



-3z F (r s ) — 5- + 
x z 

irA(r s ) cosx] + O 



[z F (r s ) since 



(26) 



see Appendix g A(r s ) = a(r s ) (^) 1 / 2 0.63 is the Fermi 
edge coefficient, the prefactor of the logarithmic term (fc— 
1) In |fc - 1| in n(k m 1, r s ). The factor 0.63 is the Kulik 
number 7.91, see Eq. dA6| ), divided by 4ir. In the inverse 
Fourier transform (0), the oscillatory terms of Eq. ( jgtj ) 
do not affect the small-A: behavior of n(k,r s ), because 
their average is zero. Since n(k <C l,r s ) = no(r s ) + 
0(k 2 ), the large- a; behavior of the non-oscillatory fix, r s ) 
is cx 1/x 6 or faster. f(x,r s ) is displayed in Fig. H. One 
may partition n(k, r s ) and correspondingly /(a;, r s ) in the 
following way: 



n(k,r s ) = z F (r s )6(l - k) + m(k,r s ), 
f(x,r.) = 3z F (r s )^ + A(x,r s ), (27) 



where j%(x) = (sin x — x cos x)/x 2 . ni(k,r s ) is a continu- 
ous function with ni(l",r s ) = ni(l + ,r s ) and an infinite 
slope at k = 1. Figure |To| shows ni(k,r s ) and fx(x,r s ) 
for r s = 5. 

The 1-matrix squared appears in the cumulant parti- 
tioning of the pair density, 

g(x,r s ) = 1- -\f(x,r s )\ 2 - h(x,r s ), x = k F r 12 . (28) 

This defines the cumulant pair density h(x,r s ), which is 
the diagonal of the cumulant 2-matrix x(l|l', 2|2'). The 
spin-resolved version of Eq. (28) is 

9n(x,r s ) = 1 - \ f(x,r s )\ 2 - h r i(x,r s ), 
9n( x >i~s) = 1 - h n (x,r s ) (29) 

with g{x,r s ) = \[g^{x,r s ) + g n (x,r s )} and h(x,r s ) = 
\[h^{x,r s ) + hfi(x, r s )]. Notice that the (generalized 




FIG. 10: Continuous part of n(k,r s ), ni(k,r 3 ) = n(k,r s ) — 
2f(j"s)0(1 — k) (upper panel), and the corresponding 1-matrix 
fi(x,r s ) (lower panel). The oscillations of fi(x) are due to 
the infinite slope of ni(k) at k = 1. 



exchange or) Fock term |/(a;,r s )| 2 appears only in the 
parallel-spin pair density and not in the antiparallel-spin 
pair density. g^{x,r s ) describes the Fermi hole (due 
to both Pauli and Coulomb repulsion) with g^(Q,r s ) = 
h^(0,r s ) — and gn(x,r s ) describes the Coulomb hole 
(only due to the Coulomb repulsion) with g(Q,r s ) < 1. 
In addition to the above mentioned correlation-strength 
indices, the quantities h"^(0,r s ) measuring the on-top 
Fermi-hole curvature, and /i||(0, r s ) measuring the on- 
top Coulomb hole are other ones. 
From Eq. (p6|) it follows 



|/(*»l,r s )| 5 



z F (r s ) 2irA(r s )z F (r s ) 



-9—^-= — sin 2x 



O 



(1 + cos 2x) 
(30) 



If this is inserted into Eq. (f28|), then the non-oscillatory 

terms, + 97rz g A , are canceled by the asymptotics of 
h(x, r s ), which follow from the sum-rule properties of the 
static structure factor S(q,r s ), see Sec. and Ref. 
The nominator of the oscillating 1 /x 5 term can be written 
as 



9z F ^J z 2 {r s ) + TT 2 A 2 (r s ) cos(2a; + 2x (r s )), 
z F (r s ) 



tan 2xo(r s ) 



irA{r s 



8 




1 
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0.6 
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0.2 




x - k F r 12 



r s = 5 



4 

x = 



k F r 12 
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FIG. 11: Cumulant pair densities for parallel and antiparal- 
lel spins, obtained by combining the present work with the 
results of Ref. |!e| . 



(from which it follows that h 2 (q,r s ) has a discontinuity 
in its second derivative at q = 2) and the cumulant com- 
ponent 

h(q, r.) = a 3 [°° dx 3 S -^h(x, r s ), (33) 
Jo 1 X 

which is simply the Fourier transform of the cumulant 
pair density h(x,r s ). In Eq. ([52]) the convolution theorem 
has been applied. n 2 (q, r s ) is related to the probability of 
finding a pair of electrons with given relative momentum 
qfflB 

Notice that the sum rule S(q — + 0, r s ) — is equivalent 
to the sum rule 



dx 3 h(x, r s ) 



dk 3 n(k,r s )[l-n(k,r s )}. 



(34) 

The l.h.s. equals h(0,r s ) and the r.h.s. equals 1 — 
^n 2 (0,r s ). It was already P--CX Lowdin who has asked 
what meaning the r.h.s. has .t3 According to Eq. (|34|), 
it fixes the normalization of the cumulant pair density 
h{x,r s ) and is another particle-hole symmetric measure 
of the correlation strength; it is also reported in Fig. ||. 
Eq. (Q) is sometimes called perfect screening sum rule 
or charge neutrality condition. 

For non- interacting electrons (r s = 0), the cumulant 
part vanishes, h(q,r s ) = 0, and the Fock part S F (q,r s ) 
simply yields 



The on-top properties of g-f\(x -C l,f s ) and 
g^l(x <c pljrPx) are determined by the coalescing cusp 
theorems J£3ti3 

With Eq. (p9|), with the spin-resolved pair densities 
of Ref. [l5], and with f(x, r s ) of this paper, the resulting 
cumulant pair densities h-\-\(x, r s ) and hn(x, r s ) are plot- 
ted in Fig. P 
those of Ref. 



For small r s (<C 1), our results agree with 



So(q) = 



<l 



1 (1 



2 V2 



1 - 



<l 



- 1 



(35) 



with the linear small-g behavior 3g/4 . . Fo r interacting 
electrons (r s ^ 0) the small-g sum ruleBo 



S{q<^ l,r s ) = 



1 



2(ar s ) 2 ui p i(r s 
and the large-g sum ruleEl 



+ 0(q% (36) 



V. STATIC STRUCTURE FACTOR, DENSITY 
FLUCTUATIONS, AND POTENTIAL ENERGY 



S(q > l,r s 



, , 1 / 1 

1 - — ar s g (r s ) • ^ + O I ^ 



3tt 



(37) 



The cumulant partitioning of Eq. (28) causes corre- 
sponding decompositions of all the quantities containing 
l—g(x, r s ). Such quantities are the static structure factor 
S(q,r s ), the fluctuation AA r o(r s ) of the particle number 
in a fragment fl, and the potential energy v(r s ). 

So, the static structure factor is given by 

S(q, r s ) = 1 - in 2 (q, r.) - h(q, r s ) (31) 
with the (generalized exchange or) Fock component 

n\q,r s ) = a 3 dx 3 S ^\f(x, r s )| 2 = 
Jo 1 X 

+1 



J dk 3 7i{k,r s ) J dCn(y/k* +q 2 - 2kqC,r s \i2) 



hold. Wp,(r s ) = 4ire 2 p/m — 3/r 3 a.u. defines the plasma 
frequency. 

The non-idempotency and the singularities of n(k, r s ) 
determine the small-g behavior of SF(q,r s ) = 1 — 
\n 2 {q,r s ), 

3 

5 F (g<l,r s ) = S F {0,r s ) + -zl{r s )q- A(r s )z F (r s ) x 



q 2 \nq + 0(q 2 



(38) 



as shown in Append ix [d[ Notice that 5f(0, r s ) is equal 
to the r.h.s. of Eq. (|34|). In Fig. |l|, we report S F (q,r s ) 
for the ideal gas (r s — 0), for r s = 5, and in the WC 
limit. 

Eq. (|§), together with the sum rule of Eq. @, al- 
lows to extract the large- a; behavior of h(x,r s ). Namely, 




, s - (ideal gas) 
r„ = 5 



1.5 

q 



FIG. 12: The Fock component of the static structure factor 
for r s = (ideal Fermi gas), r a — 5 (present work), and in 
the Wigner crystal (WC) limit for r 3 = 75. 



because of Eq. (|3q), h(q,r s ) must cancel both the linear 
and the <7 2 lng term of n 2 (q,r s ). This implies that the 
large-x behavior of the oscillation-averaged (h(q,r s )) is: 



cumulant partx 



2 roo 



v{r s ) = -- 



dx 3 



i s Jo 

a 2 3 



-\f(x,r s )\ 2 + h(x,r s ) 



r s 2 j 



dq 3 



n 2 {q,r s ) + h(q,r s ) 



1 

x 

^ 42 ) 



(in ryd). The Fock part can also be written as@ 
^F(r s ) = — - / dkn(k,r s ) / dk'n(k',r s ) x 



kk! In ■ 



(43) 



In lowest order, with n(k,r s ) — > 9(1 — k), Eq. (43) yields 
—3/(2n ar s ). The logarithmi&iterm of v(r s —* 0) arises 
from vc(r s ), not from t!F(r s ).l 



VI. SUMMARY AND OUTLOOK 



l,r s )) =-- 



9z F (r s ) 9 7rA(r s )z F (r 5 ) 



_ (39) 
In Eq. ( pq ) these terms cancel with the non-oscillatory 
long-range part of \\f(x, r s )\ 2 . 

If we consider within the uniform electron gas a certain 
fragment f2 (e.g., a sphere of radius R) containing on 
average Nq = fi/(47rrf/3) = (R/r s ) 3 electrons, and ask 
fcj£-|hj} particle-number fluctuation AAq then the answer 



(AjV n ) 2 



^fj Q d 3 Xl J n d 3 x 2 



±\f(x,r s )\ 2 +h(x,r s ) 



= 1 - 



(ar s f 3 



WfarJ+hfar,) (40) 



Again one may ask how differently the Fock and the cu- 
mulant parts contribute to their sum and-itOithe conclu- 
sion 'correlation suppresses fluctuations' Bc3 In the case 
of a sphere SI = 4irR 3 /3, the term in the modulus in 
Eq. @ is just Q3ji(qR)/(qR), so that the Fock term 
yields 



(AN, 



°- ) 



Vo(ar s ) 9 H dq 3 
* Jo 



3.]l(qR) 
qR 



n%r s )(41) 



In Rcf. [30| it was shown that the convex Kulik function 
G(x), decorated with appropriate prefactors and with an 
appropriate inhomogeneous scaling of its argument, re- 
produces the momentum distribution n(k,r s ) of the un- 
polarizcd uniform electron gas of density p = 3/47rr 3 in 
the metallic-density regime r s € [1,6]. The r s -functions 
n(0, r s ), 71(1^ r s ), the on-top pair density g(0,r s ), and 
the kinetic energy t(r s ) form the input for such construc- 
tion. In this work, we improved the parametrization of 
n(k,r s ) via the Kulik function, and we extended it up to 
r s < 12, including the high-density regime [Eqs. (p|)-(L2|), 
and Fig. || . 

The Fourier transform of n(fe,r s ) yields the one-body 
reduced density matrix f(x,r s ) (Figs. ^ and |ic| ), with 
large x-oscillations arising from the Fermi gap zp(r s ) 
and the Fermi edge coefficient A(r s ), the prefactor of 
the logarithmic term in n(k ss l,^), which is included 
in our parametrization for the first time (Fig. ^). Sev- 
eral measures of the correlation strength have been dis- 
cussed (Fig. ||). With reliable models for the pair den- 
sity, g(x,r s ), the cumulant pair density h(x, r s ) = 1 — 
l|/(x,r s )| 2 — g(x,r s ) has been extracted (Fig. [ll]), as a 
prestep of its diagonalization in terms of cumulant gemi- 
nals (analog with the diagonalization of the pair density 
in terms of Overhauser geminals). Future work also in- 
cludes the generalization to the partially polarized gas. 
In this case, with ( = (N-\ — N±)/N, one has to con- 
sider different cases. For spin polarization £ between 
and 1, two momentum distributions are to be described, 
n-[(k,r s ,() for the spin-up electrons and ni(k,r s ,Q for 
the spirtdown electrons. So far, only the input data 
go(r s , C)E3 an d t( r s, C)li3 are available in this more general 
case. 

A small FORTRAN subroutine, which numerically 
evaluates our parametrized n(fc, r s ), is available at 



Also, the potential energy v(r s ) consists of a Fock and a http://axtnt2.phys.uniromal.it/PGG/elegas.html 
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APPENDIX A: RANDOM PHASE 
APPROXIMATION 

In RPA it isllEI n(k,r s ) = 1 - {ar s / V 2 ) 2 H (k , 1) for 
k < 1 and (ar s /TT 2 ) 2 H(k 1 1) for k > 1, where 



and 



H(x < l,y) 



Q{q,u) = 2ir\ 1 + 



1 

x 



1+X dq f°° , 
— / du 

i-x 1 Jo 



l-x z 
2q 



Q(q,u) 



y^Q(q,u) 



dq 

l+x 1 Jo 



du 



H(x>l,y) = - 



x+l ^ 

— / du 



x Jx-i 1 Jo 



2q 



(§ + x) 2 +u 2 

rr 1 



Q(q, u) 



(^) 2 +u 2 (x-l) 2 + u 



q 2 + y^Q(q,u) 
q 2 + V^Q(q,u) 



l + w 2 -g 2 /4, (l + q/2) 2 +u 2 

- In r-r^ „ — u 



2q 



(1 - q/2) 2 + u 2 



arctan I±^ + arctan 1^ 



For small or large k (far from the Fermi edge) it is 
H(k,l) -> F(k) with F(k) = H(k,0). F(k) has the 
small and large k properties 



F(k « 1) = 4.112335 + 8.984373 • k 2 + 0{k A ) (Al) 



and 



F(k»l) = — -- + 0[ — 



(A2) 



respectively .013 The coefficient of 1/fc 8 is 8.77298. For 
k near the Fermi edge it isEErEj 



H(k,l 



TT 2 1 

2ar s k 2 



G 



\k-l\ 



with 

G(x) 

and 



du 



R'(u) u R(u) - R(y) 



(A3) 



R(u) u + y u — y 



R(u) = 1 — u arctan — . 



G(x) has the small- a; behaviorEH 

G(x < 1) = G(0) + [tt (f + Vs) ■ x + Q{x 2 ) 



y= x / \/R(u) 
(A4) 



(A5) 



lnx + 0(x) 
(A6) 



with 
G(0) 



du h^^Jrr ■ arctan - w 3.353337. (A7) 
R(u) u 



The coefficient of a; In a; is 7.908799 (the Kulik number). 
G(x) has the large-x behavior 



1 



G(a:»l) = -(l-ln2)-- + 0( - 



1 



(A8) 



The coefficient of 1/x 2 is 0.160668 (the Macke number). 
The Kulik function G(x) is shown in Fig. [l]. 



APPENDIX B: THE MOMENTUM 
DISTRIBUTION OF THE WIGNER CRYSTAL 

In the low-density (large_r s ) or strongly correlated 
limit, the electrons localizcO and form a ferromagnetic 
body-centered cubic lattice with an_. electrostatic (or 
Madelung) energy of — 1.792/r s rydj£3 The next term, 

+2.65/rs^ 2 ryd, describes the coupled ,-hp^iionic zero- 
To esti- 



temperature motion in lowest order 
mate the corresponding n(k, r s ), we define with 'dhw/2 = 

2.65/ri^ 2 (in ryd, or ui = 0.88/r^ 2 in a.u.) the fre- 
quency of indepent oscillating electrons (Einstein model) . 
So, from the momentum distribution of the harmonic- 
oscillator ground state it follows 



n(k, r s — * oo) 



4?r 



1 



3 (iTUj/k 2 ) 3 / 2 



(Bl) 
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see Refs. 71 73, and p. 19. Note that k is dimen- 
sionless (measured in units of kp) and that iojk\ — 

0.88 • a 2 rl /2 = 0.24 • r\ /2 . In Ref. fr| the factor 1 is 
used instead of 0.88. n(k, r s — > oo) is correctly normal- 
ized and yields with Eq. (0) the kinetic energy (in ryd) 



t{r s 



oo) = 5 



1 2.65 



3/2 



(B2) 



as it should. The corresponding potential energy is 

1.792 1 2.65 
r, + 2 



w(r s 



oo 



3/2 



The inflexion-point trajectory with r s 
described by (see left panel of Fig. 0) 



(B3) 
as parameter is 



hnz(r s ) = {u J /{2k 2 )) 1 ' 2 = 0.35 • r\l\ 
4tt e- 1 ' 2 3.88 



n in fi(r s ) 



3 (vrw/fc|)3/2 ; . 



3/4- 



(B4) 



The region k > /c; n fl (maybe to be referred to as correla- 
tion tail) contributes to the normalization the constant 
amount 



(B5) 



r°° "1/2 
/ dk 3 n{k,r s ) = Erf (— =) - \ — = 0.80 
Jk lna V2 V en 



From Eq. flBlj ) it follows for n(0, 
4 / 1 



— * oo) 

3/2 



n (r s 



oo 



tor 1 / 3 



0.24r. 



1/2 



6.40 

1574' 



(B6) 



see Fig. |[ and for the curvature at the centre (the coef- 
ficient of k 2 ) 



26.71 



3^/2 ( .88a 2 r s 1/2 )5/2 " r V4 



(B7) 



see Fig. § 

A more refined treatment takes into account that 
in harmonic approximation there are two transversal 
branches of harmonic lattice vibrations, uj tl 2 (q,r s ), and 
one longitudinal branch u>i(g,r s ) in the face-centred cu- 
bic Brillouin zone, satisfying the sum rule cof (q, r s ) + 
u> 2 2 (q, r s ) + uf(q,r s ) = w^(r 8 ). For q = it is 
w tl 2 (0,r s ) = and therefore w;(0,r s ) = u p i(r s ). But 
also in this case the virial theorem holds and ass uming 
that n(k,r s ) is a Gaufi distribution, then Eq. (Bl) turns 
out again. F or p honons in Wigncr crystals near melting 
see also Ref. [74|. 



APPENDIX C: 1-MATRIX NEAR THE 
DIAGONAL AND FAR FROM IT 

The equation -y(l]l') = p5 aia ' i /(&f|iIi — nil) defines 
the dimensionless 1-matrix f(x). Its small-a: behavior 



of Eq. ( |25| ) follows from the large-fc behavior of n(k) 

[Eq. @]. Namely, with siny/y = 1 - y 2 /3! + y 4 /5! 

and with the integrability of n{k)k 2 k v for v = 0, . . . , 4 it 
is 



dk 3 n(k)k , 

\dy 



smy 



(CI) 



y=o 



which yields the first three terms of Eq. (|25|). Here 
/ W (0) = (£)"/WI,=o. 

Because n(k)k 2 k 5 is non-integrable, onaiias to com- 
pute /^(O) with the Kimball procedure,E3E3 which de- 
fines by 



n(fc) 



(l + /c 2 ) 4 



C = ^2 ( ar ^) 2 .9o(^) 



(C2) 

a stronger (namely ~ 1/fc 10 for fc — > oo) decaying func- 
tion Af(k), so that Af(k)k 2 k 5 is now integrable, yielding 
because of (d/dy) 5 siny/y\ y =o = 0. Thus 



/ (5) (o)=m 



(C3) 



x=0 



g?/c 3 sin /ex „ 7T _ 9n 
7TTW^ = C 52 (3 + 3s + ^ e " 



/ (5) (0) 



-C- 



9tt 



(C4) 



K) 2 !?o(r s ), (C5) 



with 

p(») = c 
It follows 

q.e.d. 

The large- a; behavior ( f26| ) follows from Eq. (^J) by par- 
tial integration. Thereby the discontinuities of n(fc,r s ) 
at k sa 1 appear. They determine the amplitudes of the 
Friedel oscillations: 

3 cos 3- 3 sin x 

f(x,r s ) = -z F (r 8 ) 5 !~ 2 F(r s ) s h fi(x,r 8 ) 

X x° 

(C6) 



with 

/i(x> l,r.) 



3 7T 

-A(r a )-5-{[- + Si(af)] cosx - Ci(x) x 
x d 2 

1 



sin a;} + O ( — r 



-A(r s )7T 



x 

3 cos a; 
x 3 



O ^ 



(C7) 



zf(t s ) is the Fermi gap and A(r s ) is the Fermi edge co- 
efficient. 



APPENDIX D: FOCK COMPONENT OF THE 
STATIC STRUCTURE FACTOR AT SMALL q 

According to the definition of Eq. ( |32] ) , the oscillations 
of \ f(x > 1, r s )| 2 [see Eq. @] only affect the disconti- 
nuities in h 2 (q,r s ) and in its derivatives at q — 2, while 
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the small-g behavior of n 2 (q,r s ) is only affected by the 
oscillation-averaged part of \ f(x, r s )\ 2 , i.e., by 

(/»(.» l,r.)) = |^ + 9«M + (^ 

ssa <D1) 

Following the procedure of KimbalBLj we define a func- 
tion T(x) by 



so that T(x — > oo) oc l/x 6 . Then, the second term will 
give the coefficient of the linear term in n 2 (q,r s ), while 
the third term will give the coefficient of a term oc q 2 In q. 
By carrying out the calculations one obtains Eq. ([3g|), 
q.e.d. . 



9 z U r s) , ^A{r s )z F (r s ) 



2{l + x 2 ) 2 (1 + x) 



(D2) 
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